# Implement the rescaled change time method on single-day data
# For AAPL 2019-01-02, WITHSIZE, SUPPORT 4, DELTA 0.5

library(data.table)
library(tidyr)
library(xts)
library(dplyr)
library(pbmcapply)
#source("/home/qfwu0528/Desktop/Research/R/qianfan/hawkes_utilities.R")
#source("~/Desktop/orderbook/Research/R/gof/load_book.R")

load("~/Desktop/orderbook/Research/R/gof/AAPL_2019-01-02.RData")
support_max = 20
delta = 0.125

#book = book[1:50000,]


# Load Hawkes estimation data
hawkes_pre = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/hawkes_markov_pref_price_eventsize_",
                          data_file_name, paste0("_support",toString(14.5),"delta", toString(delta),".rds") ))
hawkes = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/hawkes_markov_pref_price_eventsize_",
                        data_file_name, paste0("_support",toString(support_max),"delta", toString(delta),
                                               "_state_hawkes_time_customizedtwo_LASSO_0005.rds") ))


delta = hawkes$delta
support_max = hawkes$support_max
support_seq = seq(0,support_max,delta)

book =  book[modifying_queue %in% -3:3]
char_cols <- c("time", "note","ref_price_event", "queue_and_event", "event_type_coarsen")
breaks_2 = c( 34200 + seq(0,1800,60),  36000 + seq(300, 19800, 300),   55800 + seq(60, 1800,60) )
book[, time_customizedtwo := cut(book$time, breaks = breaks_2, labels =c(1:126),include.lowest=TRUE,right = FALSE)]
book$time_customizedtwo = as.numeric(book$time_customizedtwo)
state <- data.matrix(book[,mget(setdiff(colnames(book), char_cols))])
event_name <- paste(colnames(event)[-1], collapse = ".")
paste0(rep(sprintf("%+d", c(1:3, -(1:3))), each=3), "(", c("+", "-", "t"), ")")
e1 <- "+1(+).+1(-).+1(t).+2(+).+2(-).+2(t).+3(+).+3(-).+3(t).-1(+).-1(-).-1(t).-2(+).-2(-).-2(t).-3(+).-3(-).-3(t)"
e2 <- "-3.-2.-1.1.2.3"
e3 <- paste(e1, "p+(+).p+(-).p+(t).p-(+).p-(-).p-(t)", sep=".")
stopifnot(event_name == e3)
time_customizedtwo <- state[, "time_customizedtwo"]
#price <- state[,"ref_price"]
#predetermined_price_unit <- event_size_by_ticker[ticker == data_ticker, price]
#if(length(predetermined_price_unit) != 1){
#  price_unit <- predetermined_price_unit
#}else{
#  price_unit <- 10000
#}
#price <- floor(price / price_unit)
state <- state[,c("q_std_plus_1", "q_std_plus_1", "q_std_plus_1",
                  "q_std_plus_2", "q_std_plus_2", "q_std_plus_2",
                  "q_std_plus_3", "q_std_plus_3", "q_std_plus_3",
                  "q_std_minus_1", "q_std_minus_1", "q_std_minus_1",
                  "q_std_minus_2", "q_std_minus_2", "q_std_minus_2",
                  "q_std_minus_3", "q_std_minus_3", "q_std_minus_3",
                  "bid_ask_spread", # for p+(+)
                  "q_std_best_ask", # for p+(-)
                  "q_std_best_ask", # for p+(t)
                  "bid_ask_spread", # for p-(+)
                  "q_std_best_bid", # for p-(-)
                  "q_std_best_bid" # for p-(t)
)]
state[state > 9] <- 9
state_list <- lapply(1:(ncol(event)-1), function(i){
  cbind(state=state[,i])
})
state_list$time_customizedtwo = time_customizedtwo
book = book[, c("time","queue_and_event", "time_customizedtwo")]
#book[, price:=state_list[[1]][,2]]

# Selected event for goodness-of-fit testing is "+1(+)"
event_names = colnames(event)[-1]
selected_event = event_names[1]
book[, paste0("state_",selected_event) :=state_list[[1]][,1]]




# get change time
change_time = c()
change_time = pbmclapply(1:dim(book)[1], function(i){
  c(change_time, seq(book$time[i], book$time[i]+support_max, delta  ))
},mc.cores = 10)
change_time = unlist(change_time)
change_time = sort(unique(change_time))
change_time = data.frame("time"= change_time)

book_lambda = merge(book, change_time, by = "time", all = TRUE)
book_lambda = unique(book_lambda, by ="time")
# Since book has the shape right before the events. Fill NA middle values with the next existing value 
# to reveal the state value when the intensity change is only triggered by excitement step function.
# Note every book_lambda is only for a specific event. In this example, "+1(+)"
book_lambda = book_lambda %>% fill("time_customizedtwo", .direction = "updown")
book_lambda = book_lambda %>% fill(paste0("state_",selected_event), .direction = "updown")
book_lambda = data.table(book_lambda)

est = hawkes$hawkes_kernel_est
est_events = est[rownames(est) %like% "design",]
est_state = est[rownames(est) %like% "state",]
est_price = est[rownames(est) %like% "time_customizedtwo",]

stopifnot(  !   all(is.na(book_lambda[,3]))   )
stopifnot(  !   all(is.na(book_lambda[,4]))   )


p_index = which(diff(book_lambda[["time_customizedtwo"]]) !=0)
df_price = data.table( index = p_index, price = book_lambda[[3]][p_index])
df_price= rbind(df_price,  data.table(index = dim(book_lambda)[1], price = book_lambda[[3]][dim(book_lambda)[1]] ))

s_index = which(diff(book_lambda[[paste0("state_", selected_event)]]) !=0)
df_state = data.table( index = s_index, state = book_lambda[[4]][s_index])
df_state = rbind(df_state,  data.table(index = dim(book_lambda)[1], state = book_lambda[[4]][dim(book_lambda)[1]] ))

#find_price_value <- function(price){
#  est_price[rownames(est_price) %like% toString(price) ,][[selected_event]] 
#}
find_price_value <- function(price){
  est_price[rownames(est_price) == paste0("time_customizedtwo",toString(price)) ,][[selected_event]] 
}
find_state_value <- function(state){
  est_state[rownames(est_state) %like% toString(state) ,][[selected_event]] 
}

df_price[, values:= mapply(  find_price_value,price )]
df_state[, values:= mapply(  find_state_value,state )]

rm(price, state,state_list, book_lambda, book_event, book_msg, book_order,event, ref_price_nondup, 
   average_event_size,estimate_hawkes, plot_hawkes, get_theta_cls, get_hawkes_kernel_bigdt,
   get_bigdt_list, get_theta_cls_bigdt_fromfile, plot_hawkes_kernel, organize_kernel_output, 
   get_theta_cls_bigdt)
gc()

change_time = change_time$time
book = book[, c(1,2)]



stimulatee = selected_event
support_vec = seq(0,support_max,delta)




support_vec_length = length(support_vec)
start = Sys.time()
print ("Start obtaining event index!")
#event_index = parallel::mclapply( 1:dim(book)[1], function(i){
#event_index = pbmclapply( 1:dim(book)[1], function(i){
event_index = pbmclapply( 1:20, function(i){
  print (i)
  vec = book$time[i]+support_vec
  end_index = match(vec[support_vec_length], change_time )
  start_index = match(vec[1], change_time[1:end_index] )
  match(vec, change_time[start_index:end_index]        ) + start_index -1
  
}, mc.cores = 8   )
Sys.time() - start

event_size = "eventsize"
saveRDS(event_index, paste0("~/Desktop/orderbook/Research/R/gof/images/event_index_",data_file_name,"_support",toString(support_max),"delta", toString(delta),"_", event_size,".rds"  ))

event_index = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/event_index_",data_file_name,"_support",toString(support_max),"delta", toString(delta),"_", event_size,".rds"  ))

print("Start adding up all intensities!")
lambda = rep(0, length(change_time))
for (i in 1:dim(book)[1]){
  print(i)
  stimulator = book$queue_and_event[i]
  est_sub = est_events[rownames(est_events)==paste0("design_matrix", stimulator),stimulatee]
  event_index_sub = event_index[[i]]
  for (j in 1:(length(event_index_sub)-1)){
    lambda[event_index_sub[j] : (event_index_sub[j+1]-1)] = lambda[event_index_sub[j] : (event_index_sub[j+1]-1)] + rev(est_sub)[j]
  }

}

print("Start obtaining price vector!")
price_vec = pbmclapply( 2:dim(df_price)[1], function(i){
  print(i)
  rep(df_price[i, values], (df_price[i,index] - df_price[(i-1),index]))
}, mc.cores = 8)
price_vec = c(rep(df_price[1,values], df_price[1,index]   ), unlist(price_vec))

print("Start obtaining state vector!")
state_vec = pbmclapply( 2:dim(df_state)[1], function(i){
  print(i)
  rep(df_state[i, values], (df_state[i,index] - df_state[(i-1),index]))
}, mc.cores = 8)
state_vec = c(rep(df_state[1,values], df_state[1,index]   ), unlist(state_vec))

lambda = lambda + price_vec + state_vec

saveRDS(lambda, file = paste0("~/Desktop//Research/R/gof/images/lambda_",data_file_name,"_support",toString(support_max),"delta",toString(delta),".rds"   )  )
